function ret=nad(data,m,k,T,p,trend,r,LE,exo,rest,st,bayes,quadr)
% p=num of levels lags
% T=num of obs in levels(size of initial data matrix)
q=repmat(0,[1,length(data(1,:)),T-1]);
   for i=2:T
q(:,:,i)=sum(data(1:i,:));
    end
s=squeeze(q(:,:,2:end));clear q
 s=reshape(s',length(data)-1,length(data(1,:)));
 s=[repmat(0,1,length(data(1,:)));data(1,:);s];
 if LE==1;
     data=s;
%  else  data=data;
  end
    if rest==1;
        data=[data(:,1:end-st-1) s(2:end,end-st:end)];
%                 data=[s(2:end,end-st:end-1) data(:,1:end-st-1) s(2:end,end)];

%     else
%     data=data;
    end
   clear s
x=[data(1:end,1:m)];
% m=num of exogenous variables
y=[data(1:end,m+1:k)];clear data
% k=number of variables
dy=diff(y);
%     dd=zeros(size(dy,1),1);
% %     dd([29-1:128,150-1:220],:)=1;
%     dd([5:57,141-1:185,222-1:276],:)=1;
% dy=(1-repmat(dd,1,size(dy,2))).*dy;
dx=diff(x);

for i=1:p
dx_{i}=dx(p+1-i:end+1-i,:);
end
 dx_=cell2mat(dx_);
 for i=1:p-1
dy_{i}=[dy(p-i:end-i,:)];
 end
 dy_=cell2mat(dy_);
if LE==1 
    dz=[ones(T-p+1,1) dx_ dy_];
%         dd=zeros(size(dz,1),1);
% % dd([29-6:128,150-6:220],:)=1;
%            dd([5:57,141-1:185,222-1:276],:)=1;
%      dz=(1-repmat(dd,1,size(dz,2))).*dz;
else
    dz=[ones(T-p,1) dx_ dy_];
end
    r_1=dy(p:end,:)-dz/(dz'*dz)*dz'*dy(p:end,:);
% r_1=first residuals
if LE==1
z_=[(1:T-p+1)' (1:T-p+1).^2' x(p:end-1,:) y(p:end-1,:)];
% dd=zeros(size(z_,1),1);
% % dd([29-6:128,150-6:220],:)=1;
%            dd([5:57,141-1:185,222-1:276],:)=1;
%      z_=(1-repmat(dd,1,size(z_,2))).*z_;      
else
z_=[(1:T-p)' x(p:end-1,:) y(p:end-1,:)];
end
% z_=dependent variable for second regressiom
if trend==0&&LE==1
z_=z_(1:end,3:end);
end
if trend==0&&LE==0
z_=z_(1:end,2:end);
end
if trend==1&&quadr==0
z_(:,2)=[];
end
if trend==1&&quadr==1
z_=z_(:,:);
end
if LE==1&&trend==1&&quadr==1
      dz=[repmat(1,T-p+1,1) (1:T-p+1)' (1:T-p+1).^2' dx_ dy_];
end 
if LE==1&&trend==1&&quadr==0
      dz=[repmat(1,T-p+1,1) (1:T-p+1)' dx_ dy_];
end 
H=eye(k-st,k); 
if rest==1
    z_=z_*H';
end
r_0=z_-dz/(dz'*dz)*dz'*z_;
% r_0=second residual
s_00=r_0'*r_0*(1/(T-p));
s_01=r_0'*r_1*(1/(T-p));
s_11=r_1'*r_1*(1/(T-p));
s_10=r_1'*r_0*(1/(T-p));
 [b,d]=eig((s_00)\s_01/(s_11)*s_10);
 d=diag(sort(diag(d),1,'descend'));
 if rest==1
for i=1:k+trend+quadr-st
max_eig(i)=-(T-p)*log(1-d(i,i));
% trend=1 if there's a trend, c=0 if there's no trend
% max_eig=statistic for null hypothesis of i-1 coint vectors
end
else
   for i=1:k+trend+quadr
     max_eig(i)=-(T-p)*log(1-d(i,i));
   end
 end
if rest==1
for i=1:k+trend+quadr-st
trace_stat(i)=sum(max_eig(i:end));
end
else
   for i=1:k+trend+quadr
           trace_stat(i)=sum(max_eig(i:end));  
   end
end
ret.b=b;
ret.trace_stat=trace_stat;
ret.max_eig=max_eig;
% trace_stat=statistic for null hypothesis of i-1 coint vectors
zero1=zeros(k);
b_=b(1:end,1:r);
if r==0 
 b_=zero1(1:k-m,1:r); 
else b_=b(1:end,1:r);
end
% r=cointegration rank, b_=beta matrix
if r<1
    b_nor=zero1(1:k-m,1:r);
else
b_nor=b_/(b_(1+trend:r+trend,1:r));
end
% if rest>0
%     b_nor=[b_nor;zeros(st,r)];
% else
%     b_nor=b_nor;
% end
% DIFFERENT NORMALIZATION
%   if r<1
%     b_nor=zero1(1:k-m,1:r);
% else
% b_nor=b_*inv(b_([1 6],1:r));
% end
    
ret.b_=b_nor;
% COINTEGRATION RESTRICTIONS (13 VARIABLES)
% if r<1
%     b_nor=zero1(1:k-m,1:r);
% else
%     b_nor=[0 0 1 zeros(1,6) -1 0 0 0;zeros(1,1) -1 0 0 1 0 0 0 1 0 0 -1 0]';
% end
% COINTEGRATION RESTRICTIONS (10 VARIABLES)
% if r<1
%     b_nor=zero1(1:k-m,1:r);
% else
%     b_nor=[1 zeros(1,5) -1 0 0 0;zeros(1,2) 1 0 0 0 0 1 -1 0]';
% end
% b_nor=normalized beta;
if r==0 
    a=zero1(1:k-m,1:r);
else a=s_10*b_nor/(b_nor'*s_00*b_nor);
end
ret.a=a;
% a=s_10*b_nor*inv(b_nor'*s_00*b_nor);
% a=adjustment coefficients, num of colums=r

if LE==1;
 p_=sparse(eye(T-p+1))-repmat(1,T-p+1,1)/(repmat(1,T-p+1,1)'*repmat(1,T-p+1,1))*repmat(1,T-p+1,1)';
else
 p_=sparse(eye(T-p))-repmat(1,T-p,1)/(repmat(1,T-p,1)'*repmat(1,T-p,1))*repmat(1,T-p,1)';
end
 if r<1
     G=(dy(p:end,:)')*p_*dz(:,2:end)/(dz(:,2:end)'*p_*dz(:,2:end));
 else
G=(dy(p:end,:)'-a*b_nor'*z_')*p_*dz(:,2:end)/(dz(:,2:end)'*p_*dz(:,2:end));
 end
 if LE>0&&trend>0
     G=(dy(p:end,:)')*p_*dz(:,2:end)/(dz(:,2:end)'*p_*dz(:,2:end));
     G=G(:,1:end);
 end
 www=length(G(:,1));eee=length(G(1,:));
 if exo==1
     G=zeros(www,eee);
 else
     G=G;
 end
     ret.G=G;
if LE>0
    cons=(dy(p:end,:)'-G*dz(:,2:end)')*ones(T-p+1,1)/(ones(T-p+1,1)'*ones(T-p+1,1));
else
 if r<1
    cons=(dy(p:end,:)'-G*dz(:,2:end)')*ones(T-p,1)/(ones(T-p,1)'*ones(T-p,1));
else
    cons=(dy(p:end,:)'-G*dz(:,2:end)'-a*b_nor'*z_')*ones(T-p,1)/(ones(T-p,1)'*ones(T-p,1));
 end
end
if exo==1
    cons=zeros(length(cons(:,1)),length(cons(1,:)));
end
ret.c=cons;
if LE==1;
    res=dy(p:end,:)'-(cons*ones(T-p+1,1)'+G*dz(:,2:end)');
%       res=dy(p:end,:)'-(cons*[(1-dd).*ones(T-p+1,1)]'+G*dz(:,2:end)');
   
else
if r<1
    res=dy(p:end,:)'-(cons*ones(T-p,1)'+G*dz(:,2:end)');
else
    res=dy(p:end,:)'-(cons*ones(T-p,1)'+G*dz(:,2:end)'+a*b_nor'*z_');
end
end
  if bayes==0   
clear dz dy p_ z_ cons p_ dx
  end
if LE==1&& m==k-1&&r==0;
cov_res=res*res'*(1/(T-(p-1)-(p-1)*k-1));
ret.cov_res=cov_res;
end
if LE==1&& m<k-1;
cov_res=res*res'*1/(T-(p-1));
ret.cov_res=cov_res;
end
if LE==0&& m==k-1&&r==0;
    cov_res=res*res'*(1/(T-p-p*k-1));
ret.cov_res=cov_res;
end
if LE==0&& m<k-1;
    cov_res=res*res'*1/(T-p);
ret.cov_res=cov_res;
end
if LE==0&& m==k-1&&r>0;
    cov_res=res*res'*1/(T-p);
ret.cov_res=cov_res;
end
% if LE<1&& m<k-1&&r>0;
%     cov_res=res*res'*1/T;
% ret.cov_res=cov_res;
% end

coint=a*b_nor';
if rest==1
    coint=[coint zeros(length(coint),st)];
end
%  ret.trace_stat=trace_stat;
%  ret.max_eig=max_eig;
 ret.d=d;
% if r<1
%     ret.ec=[];
% else
%  ret.ec=z_(:,1:k+trend)*b_nor;
% end
if LE==1&&trend==1&&quadr==1
     trend_coef=G(:,1);trend_coef1=G(:,2);
     G=G(:,3:end);
     ret.G=G;
     ret.trend_coef=trend_coef;  
     ret.trend_coef1=trend_coef1;
end
if LE==1&&trend==1&&quadr==0
     trend_coef=G(:,1);
     G=G(:,2:end);
     ret.G=G;
     ret.trend_coef=trend_coef;
end
ret.b_nor=b_nor;
 ret.coint=coint;
 ret.res=res;
 ret.LR=(T-p)*(log(det(cov_res)));
 
 % BAYESIAN ESTIMATION
if bayes==1
   %    GENERARTNG SIGMA MATRIX;
  
         ret.SIG = iwishrnd(res*res',length(res));
  %    GENERATING COEFFICENT MATRIX;
   dz=circshift(dz,[0 -1]);
U=[G cons];
U1=reshape(U',length(U(:,1))*length(U(1,:)),1);
   COEF=mvnrnd(U1,kron(ret.SIG,inv(dz'*dz)));
   COEF=reshape(COEF',length(U(1,:)),length(U(:,1)))';
   ret.G=COEF(:,1:end-1);ret.c=COEF(:,end);
   if exo==1
     G=zeros(www,eee);
         ret.c=zeros(length(cons(:,1)),length(cons(1,:)));
      ret.G=G;
   end
if LE==1;
    ret.res=dy(p:end,:)'-(ret.c*ones(T-p+1,1)'+ret.G*dz(:,1:end-1)');
else
if r<1
    ret.res=dy(p:end,:)'-(ret.c*ones(T-p,1)'+ret.G*dz(:,1:end-1)');
else
    ret.res=dy(p:end,:)'-(ret.c*ones(T-p,1)'+ret.G*dz(:,1:end-1)'+a*b_nor'*z_');
end
end
cov_res=ret.SIG;
ret.cov_res=cov_res;
else
    ret.SIG=0;
end
% US
%  ret=nad(data,0,7,373,6,0,2);
% sig_x=ret.cov_res;
%  coint_x=ret.coint;
%  E=ret.G;
% res_x=ret.res;
% ret1=nad(data,7,8,373,6,0,0);
% G=ret1.G;
%  coint_y=ret1.coint;
% sig_y=ret1.cov_res;
% res_y=ret1.res;
% US-FIRST SAMPLE
%  ret=nad(data,0,6,335,6,0,2);
% sig_x=ret.cov_res;
%  coint_x=ret.coint;
%  E=ret.G;
% res_x=ret.res;
% ret1=nad(data,6,7,335,6,0,0);
% G=ret1.G;
%  coint_y=ret1.coint;
% sig_y=ret1.cov_res;
% res_y=ret1.res;
%  UK
%  ret=nad(data,0,4,370,6,0,0);
% sig_x=ret.cov_res;
%  coint_x=ret.coint;
%  E=ret.G;
% res_x=ret.res;
% ret=nad(data,4,10,370,6,0,4);
% G=ret.G;
%  coint_y=ret.coint;
% sig_y=ret.cov_res;
% res_y=ret.res;
% CANADA
%  ret=nad(data,0,4,228,6,0,0);
% sig_x=ret.cov_res;
%  coint_x=ret.coint;
%  E=ret.G;
% res_x=ret.res;
% ret=nad(data,4,10,228,6,0,4);
% G=ret.G;
%  coint_y=ret.coint;
% sig_y=ret.cov_res;
% res_y=ret.res;
% SWEDEN
%  ret=nad(data,0,4,159,6,0,0);
% sig_x=ret.cov_res;
%  coint_x=ret.coint;
%  E=ret.G;
% res_x=ret.res;
% ret=nad(data,4,10,159,6,0,4);
% G=ret.G;
%  coint_y=ret.coint;
% sig_y=ret.cov_res;
% res_y=ret.res;